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SUMMARY 

Using  a  generalization  of  Newton's  method,  a  nonlinear 
parabolic  equation  of  the  form  u.  —  u  ■  g(u),  and  a 

v  XX 

nonlinear  elliptic  equation  u  +  u  -  eu,  are  solved 

xx  yy 

numerically.  Comparison  of  these  results  with  results 
obtained  uBlng  the  Picard  iteration  procedure  show  that  In 
many  cases  the  quasllinearlzatlon  method  offers  substantial 


advantages  in  both  time  and  accuracy. 


P-2200 

1 


SOME  NUMERICAL  EXPERIMENTS  USINO  NEWTON'S  METHOD  FOR 
NONLINEAR  PARABOLIC  AND  ELLIPTIC 
BOUNDARY-VALUE  PROBLEMS 

Richard  Bellman 
Mario  Juncosa 
Robert  Kalaba 


1 .  INTRODUCTION 

The  numerical  treatment  of  Initial— value  problems  In 
ordinary  differential  equations  on  an  electronic  digital 
computer  usually  is  no  more  involved  in  the  nonlinear  case 
than  in  the  linear  one.  In  the  handling  of  boundary— value 
problems  this  is  not  so.  In  the  linear  case,  when  a  solution 
exists,  the  applicability  of  the  superposition  principle 
provides  a  decided  advantage  in  boundary— value  problems  in 
that  it  leads  to  solving  at  most  a  number  of  Initial-value 
problems  equal  to  the  order  of  the  system.  On  the  other  hand, 
in  the  nonlinear  case,  one  of  the  possibilities  for  solution 
is  to  resort  to  an  iterative  technique  which  replaces  the 
problem  with  a  sequence  of  linear  problems  in  which  one  can 
use  the  superposition  principle. 

In  solving  the  differential  equation 

(1)  L(u)  -  f(u) 

where  L  is  a  linear  ordinary  differential  operator  of  at 
least  second  order  and  conditions  are  prescribed  to  be 
satisfied  by  u(x)  at  at  least  two  points,  ne  may  linearize 
by  using  Picard's  method  which  introduces  a  sequence  of 
functions  (uv  '(x)}  which  satisfy  the  same  boundary  conditions 
as  u(x)  and  the  linear  ordinary  Inhomogeneous  differential 
equation 
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(2)  Lu(k+1)  -  f(u(k)). 


When  the  eequenoe  (u^k^(x))  converges  Ohe  convergence  is 
linear,  i.e. 


(3)  u(k+1)  -  u  -  0(u^k^  -  u) 


ss  k  -»  oo  . 

However,  if  f(u)  is  differentiable,  we  can  linearize 
in  a  different  way.  If  one  replaces  the  right  hand  side  of 
(2)  by 

f(u(k))  ♦  (u(k+1>  -  u^k^  )f  •  (Jk)) 

then  we  have  again  a  linear  inhomogeneous  ordinary 
differential  equation 

(4)  Lu-  f '(u(k))u(k+1)  -  f(u(k))  -  u^f '(u^) 

which  results  in  a  sequence  (u^k^(x)}  which,  when  convergent, 
is  usually  quadrat ically  convergent,  i.e. 

(5)  u^*1)  -  u  «  o((u^k^  -  u)2} 
as  k  -*  oo  . 

The  idea  of  the  use  of  second— order  convergent  iterative 
procedures  for  solving  systems  of  equations  other  than  alge¬ 
braic  or  transcendental  is  not  new.  However,  except  for  the 
work  of  Hestenes  [1]  and  Stein  (2 j  and  a  brief  mention  in 
Kline's  book,  sec.  49  [}],  it  seems  to  have  received  only 
soant  attention  in  tne  American  literature  on  numerical 
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analysis  compared  to  that  In  foreign  publications.  In  1905* 
only  a  few  decades  after  Picard's  work,  Chaplygin  [4]  presented 
what  amounts  to  Newton's  method  for  approximate  integration  of 
differential  equations.  More  recently  general  functional- 
analytic  treatments  of  Newton's  method  and  some  of  its  variants 
in  a  Banach  space  have  been  given  by  a  number  of  authors, 
notably  Kantorovich  [5*6],  Zagadskii,  Mysovskikh,  Fenyo, 

Collatz,  Schroder  [7],  Bartle,  and  Stein. 

References  to  most  of  these  and  to  others  can  be  found  in 
[2,6,7]. 

From  the  practical  point  of  view,  in  spite  of  all  this 
analytic  treatment,  because  of  the  effects  of  truncation 
errors,  round-off  errors,  crude  bounds  on  the  higher  order 
derivatives,  the  labor  of  computing  higher  order  derivatives 
or  differences  for  the  nlgher  order  methods,  the  accuracy 
desired  in  the  final  approximation  to  the  solution,  and  the 
word  length  and  computing  mode  of  the  machine  available  for 
the  computation  thf  efficiency  of  the  higher  order  convergent 
methods  vs.  that  of  lower  order  methods  cannot  be  decided 
solely  on  the  basis  of  the  order  of  convergence.  An  analysis 
of  the  interaction  of  such  effects  is  generally  far  more 
difficult  than  that  of  simply  determining  orders  of  convergence. 
This  difficulty,  if  not  Impossibility,  Justifies  some  numerical 
experimentation.  Although  Kantorovich,  Collatz,  and  Schroder, 
give  some  examples  of  applications  to  eigenvalue  problems. 
Integral  equations,  and  differential  equations  (mainly  ordinary) 
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and  He stents  and  Stein  give  some  applications  to  the  calculus 
of  variations ,  the  body  of  experimental  results  In  the 
American  literature  Is  still  quite  small. 

For  ordinary  differential  equations  there  is  quite  a 
large  number  of  stable  methods  of  numerical  Integration  whose 
truncation  errors  are  of  fairly  high  order  In  the  integration 
step  else.  Consequently,  by  the  choice  of  one  of  such  methods 
the  interaction  between  the  truncation  error*,  and  the  rate  of 
convergence  of  the  numerical  solutions  of  either  ^2)  or  of  (4) 
as  s  function  of  k  can  be  kept  quite  small.  Hence  the 
computational  effort  to  obtain  numerical  approximations  to  the 
solution  of  (1)  is  largely  determined  by  whether  ona  chooses 
to  solve  the  linear  differential  equations  (2)  or  (4),  thus 
comparing  quadratic  convergence  (5)  and  some  extra  computation 
In  the  evaluation  of  in  (4)  with  linear  convergence 

(3)  and  no  evaluation  of  derivatives  of  f  in  (2).  This, 
indeed,  has  been  oompared  by  one  of  the  authors  for  a  simple 
nonlinear  second  order  ordinary  differential  equation  with 
two  point  boundary  conditions.  See  [8]  wherein  is  given  a 
novel  relatively  general  derivation  of  (4)  utilising  the 
operatien  of  function  maximisation.  See  also  [12]  for  more 
general  applications. 

For  elliptic  and  parabolic  partial  differential  equations, 
however,  the  stable  numerical  methods  commonly  used  for 
solution  usually  have  truncation  errors  of  low  order  in  the 
mesh  sizes.  Furthermore,  the  numerical  treatment  of  elliptic 
and  multl-epace-dlmenaional  parabolic  cases  almost  always 
result  in  systems  of  algebraic  equations  which  are  the 
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linearized  discrete  analogs  of  the  partial  differential 
equations  and  are  usually  solved  by  an  Iterative  procedure 
such  as  relaxation  or  over  relaxation.  Consequently,  regard¬ 
less  of  whether  one  chooses  (2)  or  (4)  for  the  linearization 
of  (l),  there  Is  considerable  Interaction  between  the 
truncation  errors  of  the  discrete  analog,  the  rate  of 
convergence  (3)  or  (5)»  depending  on  the  choice  of  either  (2) 
or  (4),  and  the  rate  of  convergence  of  the  relaxation  or 
overrelaxation  procedure. 

As  «.  result  an  analytic  representation  of  realistic 
bounds  on  the  total  error  Is  very  difficult,  if  not  Impossible, 
to  achieve .  The  purpose  of  this  note  is  to  present  the 
results  of  some  experiments  using  (2)  and  (4)  on  two  cases 
each  oi  a  parabolic  partial  differential  equation  and  an 
elliptic  one.  These  cases  In  each  type  differing  only  in  the 
boundary  conditions  clearly  show  how  markedly  the  superiority 
of  the  method  (4)  over  (2)  is  affected  by  the  above-mentioned 
Interaction  as  the  boundary  conditions  are  changed. 

2.  THE  PARABOLIC  CASE 

In  the  Interest  of  simplicity  the  experiments  on  the 
parabolic  case  were  carried  out  on  the  numerical  solution  of 

(6)  Lu  -  ut  -  u^  -  (1  +  u2)(l  -  2u) 

over  two  different  triangles i  0  £  t  £  1  —  x,  0£x£l,  and 
0  £  t  £  1.5  -  x,  0  £  x  £  1.5*  The  boundary  conditions  In 
each  case  were  so  chosen  as  to  give  the  solution 
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(7) 


u(x,t)  •  tan(x  +  t) . 


That  this  solution  Is  unique  can  be  established  easily  by 
classical  techniques . 

Solutions  of  the  differential  equation  analogs  of  the 
Picard  iteration  (2)  and  of  the  Newton  procedure  (4)  for  the 
differential  equation  (6)  were  compared.  The  numerical 
solutions  in  each  case  were  obtained  over  the  points 
(roAt,nAt)  of  a  grid  superimposed  on  the  respective  triangles. 
In  order  to  avoid  the  possibility  of  numerical  instability 
(since  we  used  Ax  »  At  •  .01  in  the  experiments)  we  used  the 
Crank-Nlcolson  difference  operator  [9 >10]  for  the  discrete 
analog  of  (6),  Thus,  if  u^  n  designates  the  solution  of  a 
difference  equation  analog  of  (6)  at  the  lattice  point 
(mAt,nAt)  in  the  triangle  and  since  an  iterative  process  is 

needed  to  resolve  the  nonlinearity  of  (6),  Lu  is  replaced  by 
u(k+l)  _  u  ^  r 

™jn+l  Vn _ 1.  Li™)  _  ?u(k+1)  +  u(k+1) 

T  L  m,n+l  +  Vl.n- 


2(A x)‘ 


,n+l 


(knJ  (kn^  (kn^ 

+  um— l,n  2um,n  +  V*-l, 


where 


denotes  the  final  number  of  iterations  to  obtain  an 


acceptable  approximation  to 


u_  „  &t  the  grid  points  on  the 
m,n 


line  t  •  nAt.  In  the  iterative  formula  two  possibilities  for 

the  function  to  replace  f(u^)  on  the  right  hand  aides  of 

(kj 

(2)  and  (4)  were  consideredj  one  was  simply  fCi^  n  )  &nd  the 
other  was  the  average  (ffi^  ♦  f(um  n  ))/2*  However, 

sinoe  it  developed  early  in  the  experiments  that  the  latter 
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representation  was  a  bit  better  than  the  former  and  since  we 
were  primarily  concerned  with  a  comparison  of  Picard's  method 
with  that  of  Newton,  we  continued  with  the  latter  only  for  the 
remainder  of  the  work,  which  is  reported  here. 

The  significant  observation  was  that  while  the  Picard 
procedure  and  the  Newton  procedure  required  about  the  same 
amount  of  work  in  the  case  of  the  triangle  0  £  t  £  1  —  x, 

0  £  x  £  1,  about  nine  times  the  number  of  iterations  required 
for  Newton's  method  were  needed  to  obtain  comparable  accuracy 
by  the  Picard  procedure  when  the  region  of  Interest  was  the 
larger  triangle  0£t£l.5  —  x,  0  £  x  £  1.5* 

The  criterion  for  acceptance  of  an  approximation  was  a 
commonly  used  one,  viz.,  that  the  maximum  relative  change  per 
line  be  less  than  a  prescribed  amount  before  passing  to  the 
next  line.  Thuo  our  requirement  was  that  •  min  k  such 
that 


max 

m 


(*)  _  u(*-l) 

m,n  m,n 

'  m  — 


1  10 


-6 


Since  in  this  problem  the  true  solution  is  known  (7),  we  could 
have  a  criterion  based  on  the  true  relative  error.  However, 
this  is  impossible  when  the  solution  is  not  known  and  we  wished 
to  simulate  such  a  condition.  A  check  was  made  on  what  the  true 
relative  errors  were.  In  the  case  where  the  base  of  the 
triangle  was  1.0  the  maximum  true  relative  errors  on  each  line 
using  Newton's  method  and  using  Picard's  were  in  agreement  with 
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each  other  to  about  two  significant  figures.  These  maximum 
time  relative  errors  never  exceeded  2-9  x  10  ^  and  usually 
lay  near  the  line  t  ■  1  -  Jx/2 .  When  the  triangle  had  a  bass 
equal  to  1.5  the  solution  using  Newton's  method  was  run  only 
up  to  the  line  t  ■  0.22  and  with  Picard's  up  to  t  -  0.5* 
Again  for  the  values  obtained  there  was  very  close  agreement 
between  the  location  and  the  value  of  the  true  maximum  relative 
errors.  They  were  generally  located  near  the  line  x  +  t  -  1.4 
and  generally  did  not  exceed  1.8  x  10"^.  The  substantial 
difference  between  true  relative  errors  and  the  relative 
changes  between  successive  iterations  should  be  taken  as  a 
caution  to  numerical  analysts  to  set  bounds  in  stopping 
criteria  based  on  relative  changes  between  successive  iterations 
much  lower  than  the  desired  maximum  true  relative  errors. 

The  following  table  gives  the  comparative  numbers  of 
iterations  to  achieve  our  criterion  for  acceptance  of  an 
approximation  before  going  on  to  the  next  t-llne. 


BASE  OF 

NUMBER  OF 

TR IAN OLE 

METHOD 

ITERATIONS 

t- INTERVAL 

1.0 

Newton 

1 

! 

2 

(0.01,0.98) 

1 

(0^8,0.99) 

1.0 

Picard 

4 

(0.01,0.68) 

2 

(0.68,0.93) 

2 

(0-92. 0.99) 

__  1.5  .  j 

Newton 

2 

(0.01,0.23) 

1-5  . 

Picard 

28 

(0.01,0.50) 

TABLE  1 
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3>  THE  ELLIPTIC  CASE 

The  experiments  on  the  elliptic  case  were  carried  out  on 
the  numerical  solutions  of 


(8) 


Lu  ■ 


u  + 
XX 


u 


yy 


-  e 


u 


in  the  region  0  <  x  <  1/2,  0  <  y  <  1/4,  for  two  sets  of 
bound xry  conditions  u  •  0  in  one  case  and  u  •  10  in  the 
other.  The  equation  (8)  is  of  considerable  interest  in  some 
physical  problems  and  the  existence  and  uniqueness  of  the 
solutions  to  the  boundary-value  problems  we  have  here  are 
assured  by  the  classical  theory. 

In  each  of  the  two  cases  we  compared  the  two  methods  for 
linearization,  Picard's  and  Newton's,  which  resulted  in 
comparing  the  numerical  solutions  obtained  for 


(9) 


u 


(k+1)  +  u(k+l) 


XX 


yy 


u 


00 


and 

(10) 


u 


(k+1)  +  u(k+l) 


XX 


u 


(it) 


yy 


—  e 


u 


(k+l) 


00 


(1  - 


(k)  v 
uv  ' ) , 


respectively.  As  is  customary  in  discretizing  the  problem  for 

a  numerical  solution,  the  Laplaclan,  u  +  u  ,  in  (9)  and 

xx  yy 

(10),  was  replaced  at  interior  meshpoints  of  the  region  of 
interest  by  the  expression 


(11) 

where 


u 


m ,  n 


1b 


( u  ,  +  u  .  ♦  u  ,  +  u  , 

v  m+l.n  m,n+l  m— l,n  m,n— 1 

the  discrete  analog  of  u(mAx,nAt) 


4u 


m ,  n 


)/h 


3 

c. 


and 
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h  ■  Ax  ■  Ay,  which  In  the  experiments  was  set  equal  to  1/64. 

If  we  order  the  points  (m,n)  such  that  (ni,n)  precedes 
(m'jn')  if  n  <  n'  or  if  n  ■  n’  and  m  <  m'  and  denote 
the  resulting  set  of  15*31  numbers  u'  '  which  are  the  approxl- 

Ul  p  » 1 

mations  to  the  solution  at  (mAx,nAy)  at  the  k-th  Picard  or 

(k ) 

Newton  iteration  by  the  same  notation  as  above,  uv  ' ,  then, 
after  division  of  the  component  equations  by  the  appropriate 
factors,  via.,  the  negative  of  the  coefficients  of  the  central 
term,  n,  in  the  equations,  we  obtain  systems  of  linear 
algebraic  equations 

(12)  (I  -  I*  -  Uk)u(k+1)  -  b(k) 

to  solve.  In  (12),  I,  L^,  and  Uk  are  465  x  465  matrices,  I 
being  the  identity,  while  and  Uk  are  respectively  lower 
and  upper  triangular  matrices  appropriate  to  the  particular 
method  of  iteration,  Picard  or  Newton.  Thus,  for  our  problem, 
they  sure  transposes  of  each  other  and  have  only  two  diagonal 
lines  of  nonzero  elements;  in  the  Picard  iteration  these  ele¬ 
ments  sure  Identically  equal  to  1/4,  while  in  the  Newton  pro- 

cedure  they  are  equal  to  (4  +  h  exp  uv  ')  where  the 

( k ) 

approximation  uv  '  is  evaluated  at  the  central  point  (m,n) 
of  the  star  of  points  indicated  in  (11).  The  components  of  the 
vector  b^k^  In  (12)  are  equal  to  -hc/4  times  the  appropriate 
values  of  the  right  hand  sides  of  (9)  and  (10)  respectively. 

Since  in  many  realistic  problems  the  size  of  the  problems 
is  so  huge  as  to  preclude  the  use  of  Qaussian  elimination  to 
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solve  the  algebraic  systems  (12),  we  chose  the  successive 
overrelaxation  method  of  Young  [11],  primarily  because  it  Is 
the  simplest  iterative  method  which  is  substantially  better 
than  the  Gauss— Seidel  method  although  it  is  recognized  that 
faster  converging  block  relaxation  and  alternating  direction 
methods  are  not  much  more  complicated  than  successive  over- 
relaxation  . 

Applying  the  successive  overrelaxation  method  directly 
to  the  systems  (12)  for  a  fixed  value  of  k  would  yield  the 
iterative  formula 

(15)  (I  -  ojLj^u^1)  -  [csU^  -  (co  -  l)l]uj,k+l)  +  xb(ic), 

X'  •  0,1,2,..., 

where  cd  is  the  overrelaxation  parameter  and  the  component 
equations  are  solved  successively  in  the  order  which  produces 
the  components  of  the  vector  u£k^^  consecutively  in  the 
order  Indicated  above.  However,  It  is  clear  that  not  much 
effort  should  be  expended  in  obtaining  a  vei'y  accurate  solu¬ 
tion  to  (12)  if  u(lc+l)  is  not  a  good  approximation  to  the 
solution  of  (8).  Similarly,  there  la  not  much  point  to 
Iterating  on  the  index  k  if  the  approximations  given  by  jl}) 
are  too  rough  as  a  consequence  of  terminating  the  iteration 
on  r  too  soon.  Thus  there  exists  the  open  problem  as  to 
when  one  should  Iterate  on  k  or  on  r  at  each  step.  We 
avoided  the  attempt  at  this  relatively  difficult  analysis  and 
simply  iterated  on  each  simultaneously  using  the  formula 
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(14) 


(I  -  toL^Ju 


-  [coU.  -  (a)  -  l)l]u^  +  cob^ 


the  component  equations  being  solved  in  the  same  order 
indicated  above .  Thus  we  see  a  strong  interweave  between  the 
numerical  method  used  for  solving  either  (2)  or  (4)  and  the 
rate  of  convergence  of  the  actual  solutions  of  (2)  or  of  (4), 
the  difficulty  of  whose  analysis  dictates  experimentation. 

We  also  note  that  this  blend  of  iterative  procedures  for  the 
solution  of  the  linear  systems  (12)  and  the  Newton  or  Picard 
iterations  to  solve  the  original  nonlinear  problem  has  the 
advantage  that  one  does  not  have  to  go  to  the  full-scale 
effort  of  solving  accuretely  the  system  (12)  for  each  value 
of  k. 

As  in  the  parabolic  case  the  criterion  for  stopping  the 
iterations  was  that  the  maximum  of  the  absolute  value  of  the 
relative  change  between  consecutive  iterations  of  the  functional 
values  be  no  greater  than  lO-^.  If  thiB  were  truly  a  bound 
on  the  relative  error  of  the  numerical  solution  and  the  solu¬ 
tion  of  nonlinear  equation  (8)  with  the  Laplaclan  replaced  by 
the  expression  (ll(,  then  this  criterion  would  be  about  two 
orders  of  magnitude  too  stringent  to  be  Justified  by  the 
truncation  errors.  However,  some  measure  of  stringency  is 
dictated  by  the  heuristic  character  of  the  stopping  criterion 
when  one  knows  neither  the  true  solution  nor  effective  bounds 
on  it . 


In  the  case  of  Picard  iteration  where  and  are 


independent  of  k,  if  one  were  to  ignore  the  fact  that 


,(*) 
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depends  on  k  the  theoretically  optimal  value  of  cjd  is 
1.752  for  the  fastest  convergence.  For  the  Newton  Iteration 
the  value  would  be  a  little  smaller,  again  ignoring  the 

( ir  \ 

change  in  and  in  as  well  as  in  bv  '  from  iteration 

to  iteration,  which  we  clearly  cannot  do.  However,  in  the 
case  of  zero  boundary  conditions  for  the  Picard  iteration, 
where  and  are  constant,  the  best  value  (or  values) 

of  co  are  undoubtedly  near  the  theoretical  value  noted  above 
since  the  change  in  u'  '  over  the  region  is  very  small  (The 
value  at  the  center  of  the  region  obtained  for  the  final 
solution  is  0.00(0 fl  to  four  significant  figures.).  Conse¬ 
quently,  for  this  case  we  experimented  only  with  the  values 
co  «  1.70  and  1.74.  No  such  confidence  existed  for  the  case  of 

Picard  iteration  on  the  situation  where  u  is  equal  to  10  on 

f  k ) 

the  boundary  because  of  the  large  variations  in  uv  '  over 
the  region  which  would  undoubtedly  Introduce  large  changes  at 
a  point  from  iteration  to  iteration.  Of  course,  in  the  Newton 
iterations,  since  and  change  from  iteration  to 

iteration,  one  has  even  less  of  an  idea  as  to  the  optimal 
value  of  co.  Consequently,  in  the  remaining  cases  a  niuiter 
of  runs  for  different  values  of  uj  were  made  to  estimate  an 


optimal  value  of  oj.  Table  ^  below  summarizes  the  results. 
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BOUNDARY 

CONDITIONS 

METHOD 

VALUE 

OF  o> 

NUMBER  OF 
ITERATIONS 

VALUE  OF  u(k> 
AT  (1/4.  1/8) 

u  ■  0 

II 

Picard 

1! 

1.74 

1.70 

64 

77 

-0 .007071 

II 

II 

It 

Newton 

1.95 

238 

It 

II 

1.74 

64 

ft 

It 

It 

1.70 

77 

II 

If 

It 

1.50 

134 

It 

1.00 

426 

II 

u  •  10 

If 

Picard 

It 

1.55 

63 

5.63995 

It 

1 . 50 

53 

It 

n 

it 

1.40 

S  ~m  r- 

66 

5.63994 

it 

ti 

It 

1.30 

71 

76 

5.63993 

II 

it 

M 

1.20 

78 

5.63992 

it 

It 

1.10 

106 

5.63397 

1.00 

134 

5.63998 

it 

m 

n 

Newton 

II 

II 

1.74 

1.70 

58 

51 

5.63995 

it 

1.63 

46 

it 

M 

1.60 

42 

it 

1 1 

1.35 

41 

ii 

II 

1.50 

50 

ii 

II 

1.40 

66 

5.65996 

1.00 

148 

5.65998 

TABLE  2 
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Prom  the  table  we  observe  that  for  the  case  of  zero 
boundary  conditions  no  advantage  was  obtained  from  the  use  of 
the  quadratlcally  converging  Newton  procedure  over  the  results 
obtained  from  the  first  order  converging  Picard  procedure.  In 
fact,  the  same  values  of  gave  the  rapidest  convergence  in 

both  cases.  On  the  other  hand,  when  the  boundary  conditions 
were  raised  to  10,  implying  rapid  changes  in  u  near  the 
boundary,  the  Newton  procedure  was  more  advantageous, 
requiring  some  41  iterations  as  compared  to  55  for  Picard's. 

We  note  further  that  curiously  the  optimal  value  of  co  was 
slightly  lower  for  Picard's  procedure. 

As  an  additional  note  we  Included  the  case  of  Qauss-Seldel 
relaxation  given  by  cu  -  1  for  comparison  with  successive 
overrelaxation.  As  expected,  the  successive  overrelaxation 
method  ranged  from  two  and  one-half  to  almost  seven  times 
faster  than  Oauss— Seidel  relaxation. 

As  an  indication  of  the  order  of  time  for  a  computation, 
it  was  observed  that  on  the  RAND  Johnnlac,  a  Prlnceton-type 
machine  on  which  the  computations  were  carried  out,  it  took 
about  eleven  seconds  for  a  Newtonian  iteration. 

Another  incidental  observation  was  that  the  points  at 

( k ) 

which  the  maximum  relative  change  In  u'  '  took  place  were 
invariably  very  close  to  the  origin.  In  most  cases  it  was  at 
the  point  (1/64,  1/64). 

4_.  CONCLUSION 

We  have  observed  in  these  examples  that  when  the  sol  tion 
of  a  nonlinear  problem  has  no  steep  gradients  there  seems  to  be 


P-2200 

16 


no  particular  advantage  in  the  use  of  Newtonian  methods  over 
those  of  Picard.  On  the  other  hand,  when  steep  gradients  occur 
then  there  le  some  advantage,  greater  in  our  experiments  for 
the  parabolic  case  than  for  the  elliptic  case,  where  there  is 
considerable  Interaction  between  the  numerical  method  of 
solution  of  the  linearized  problem  and  the  convergence  rate 
of  the  Iterated  solutions  of  the  linearized  problem. 

We  remark  finally  that  it  is  possible  for  the  operator 
L  to  be  quaslllnear  as  well,  in  which  case  there  are  obvious 
modifications  to  linearize  the  nonlinear  part  of  the  operator. 

We  thank  Alfred  B.  Nelson  who  programmed  the  experiments 
for  the  Johnniac  and  Ardls  McCarroll  who  collocated  the  results . 
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